#include "PO_MoM.h"
#include "RadiationPattern.h"
#include <cmath>
#include <fstream>
#include <iostream>
void Output()
{
	std::ifstream fin("request.dat");
	double ThetaStart, ThetaEnd, PhiStart, PhiEnd;
	int ThetaNum, PhiNum, PolarType;
	fin>>ThetaStart>>ThetaEnd>>ThetaNum>>PhiStart>>PhiEnd>>PhiNum>>PolarType;
	double ThetaStep = ((1 == ThetaNum) ? 0.0:(ThetaEnd-ThetaStart)/(ThetaNum-1));
	double PhiStep = ((1 == PhiNum) ? 0.0:(PhiEnd-PhiStart)/(PhiNum-1));
	blitz::Array<double,2> FarField(ThetaNum,PhiNum);
	double max = 0;
	for (int theta=0; theta<ThetaNum; theta++)
	{
		for (int phi=0; phi<PhiNum; phi++)
		{
			FarField(theta,phi) = Farfld((ThetaStart+theta*ThetaStep)/180.0*PI,(PhiStart+phi*PhiStep)/180.0*PI,PolarType);
			if (FarField(theta,phi)>max) max = FarField(theta,phi);
		}
	}
// 	FarField = FarField / max;
	std::ofstream fout("FarField.dat");
	for (int theta=0; theta<ThetaNum; theta++)
	{
		for (int phi=0; phi<PhiNum; phi++)
		{
 			fout.setf(std::ios_base::scientific,std::ios_base::floatfield);
 			fout.setf(std::ios_base::showpoint);
 			fout.setf(std::ios_base::right,std::ios_base::adjustfield);
			fout<<FarField(theta,phi)<<std::endl;
		}
	}
	fout.close();
}
double Farfld(double theta, double phi, int polartype)
{
	using namespace blitz;
	double thesin = sin(theta);
	double thecos = cos(theta);
	double phisin = sin(phi);
	double phicos = cos(phi);
	Position3D vRk(WaveNumber*thesin*phicos,WaveNumber*thesin*phisin,WaveNumber*thecos);
	Position3D vThe(thecos*phicos,thecos*phisin,-thesin);
	Position3D vPhi(-phisin,phicos,0);
	Complex3D PO_Far(0.0,0.0,0.0);
	Complex3D MoM_Far(0.0,0.0,0.0);
	Complex3D Total_Far(0.0,0.0,0.0);
	for (int i=0; i<iEdge; i++)
	{
		int node1 = Edge(i,0);
		int node2 = Edge(i,1);
		Complex3D tmp(0.0,0.0,0.0);
		for (int LL=0; LL<2; LL++)
		{
			int pat = Edge(i,LL+2);
			int node3 = Patch(pat,0)+Patch(pat,1)+Patch(pat,2)-node1-node2;
			double signal(1.0-2*LL);

			Position3D vRho(vPatchCtr(pat)-vNode(node3));
			tmp +=vRho * signal * exp(dComplex(0.0,1.0)*dot(vRk,vPatchCtr(pat)));
		}
		PO_Far +=  Edge_length(i) * J(i) * tmp;
	}
	PO_Far = PO_Far * 0.5;
	for (int i=0; i<iWire; i++)
	{
		MoM_Far += exp(dComplex(0.0,1.0) * dot(vRk,vWire(i))) * I(i);
	}
	MoM_Far = MoM_Far * DeltaLength * UnitVector;
	Total_Far = -dComplex(0.0,1.0) *Omega*Mu/(4*PI) * (PO_Far+MoM_Far);
	if (1 == polartype)	return abs(dot(Total_Far,vThe));
	else if (-1 == polartype) return abs(dot(Total_Far,vPhi));
	else return abs(sqrt(dot(Total_Far,Total_Far)));
}
